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ABSTRACT 

The amount of current carried by an electrical discharge in its early stages of growth is strongly dependent 
on its geometric shape. Discharges with a large number of branches, each funnelling current to a common stem, 
tend to carry more current than those with fewer branches. The fractal character of typical discharges has been 
simulated using stochastic models based on solutions of the Laplace equation [1, 2]. Extension of these models 
requires the use of electron distribution functions to describe the behavior of electrons in the undisturbed medium 
ahead of the discharge. These electrons, interacting with the electric field, determine the propagation of branches 
in the discharge and the way in which further branching occurs. 

This paper reports on the first phase in the extension of the referenced models, the calculation of simple 
electron distribution functions in an air/electric field medium. Two techniques are investigated. The first is the 
solution of the Boltzmann equation in homogeneous, steady state environments. The second is the use of Monte 
Carlo simulations. Distribution functions calculated from both techniques are illustrated. Advantages and 
disadvantages of each technique are discussed. 


1.0 INTRODUCTION 

An electrical discharge is often said to have a "fractal" character. This usually means that the discharge is 
heavily branched, often to the point of seeming diffuse. The branching character is important, in that it allows 
the discharge to generate more current than is possible for a single channel. In a loose sense, the tip of each 
branch may be thought of as producing a fixed amount of current, which is then funnelled back along the branch 
to the base of the discharge. 

Models have been developed which predict the branching behavior of electrical discharges [1, 2]. These 
models embody rules which determine how a discharge will grow incrementally. The rules are generally 
stochastic in character and based on static solutions of Laplace's equation. A more physical way of defining the 
growth rules is to use the local electron distribution function at the ends of the discharge branches. The 
distribution function contains information about the electrons' velocity behavior in the electric field environment 
and can be used to probabilistically define the direction of a branch's subsequent growth. It can also be used to 
predict the likelihood of a branch bifurcating into more branches. 

This paper reports on the first phase of a multi-phase effort, the methodology needed for the calculation of 
electron distribution functions in an electric field environment. Two techniques are investigated: solution of 
Boltzmann’s equation and a Monte Carlo procedure. Subsequent efforts in this area are intended to apply the 
derived distribution functions to discharge growth models and then to apply the models to discharges occurring on 
actual systems. 

In its most general form the electron distribution function, F, describes exactly the positions and velocities 
of all of the electrons in a system as a function of time, This form is impractical to deal with, however, so 
simpler, reduced distributions, usually designated f, are defined. The reduced distribution of most interest here is 
the single particle distribution function, which is derived by integrating F over the position and velocity 
coordinates of all but one of the electrons in the system under consideration. This leaves the function f(x,v,t), 
where x represents three spatial coordinates and v three velocity coordinates. 


* This work was sponsored by Harry Diamond Laboratories under Contract DAAL02-89-C-0075. 
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Various physically measurable quantities can be derived from f by taking moments over v. One of these is 
the density variation in space, given by. 


n(x,t) = noj f(x,v,t) dv 


( 1 ) 


IN 

Here no is defined to be the average density of the system y , where N is the total number of electrons and V is 

the total volume of the system. The mean velocity of the electrons (e.g., a drift velocity of the electrons in an 
electric field) is found from the first moment. 


V(x,t) 


/ v f(x,v,t) dv 
J f(x ,v ,t) dv 


( 2 ) 


The denominator of Equation 2 represents a normalization of the distribution function. In many cases f is scaled 
so that this denominator is unity. The second moment of the distribution function can be identified with a 
pressure tensor. An infinite number of moments can be defined, but they begin to lose simple physical 
meanings when one goes beyond the second. 

The physical meaning of f itself is as follows. The quantity n 0 f(x,v) dx dv is the total number of 
electrons in the volume dx dv centered on (x,v). Therefore f(x,v) is in some sense the probability of finding a 
particular electron in dx dv centered on (x,v). 

2.0 SOLUTION OF THE BOLTZMAN EQUATIONS 

The evolution of the electron distribution function is described by the Liouville theorem [3]. Although 
somewhat simplified, this theorem essentially states that, in the absence of interactions between particles, the 
total time derivative of f is zero. 


tf 

<k 


= 0 


Expanding, one finds. 
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Here a is the acceleration experienced by the electrons from external forces (e.g., electric fields, magnetic fields, 
gravity). Equation S is commonly known as the collisionless Boltzmann equation, or the Vlasov equation. 
When collisions, either with other electrons or background particles, are taken into account, a source term must 
be included on the right hand side of Equation S. The form of this source term can be quite complicated and is 

often written formally as With the addition of this source term. Equation 5 becomes the full Boltzmann 

equation. 



ar ar _ar\ 

3x + a * dv ~ dt/coll 


( 6 ) 


As was done with the distribution function itsdf, moments can be taken of the Boltzmann equation. This 
is done by multiplying both sides of the equation by v n , where n represents the desired moment, and then 
integrating over all velocity space. The zeroth order moment of the Boltzmann equation reduces to a statement of 
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conservation of particles. The first order moment reduces to conservation of momentum, and the second order 
moment to conservation of energy. 

A form for the collision source term in Equation 6 can be derived if one makes some assumptions about the 
physical situation. For example, in a fully ionized gas, collisions are generally coulombic in nature, and 
therefore long range. Under these conditions Fokker-Planck theory is appropriate for calculating transport 
properties. Alternatively, the case considered here is one of a weakly ionized plasma in which nearly all 
interactions take place between electrons and neutrals. These collisions can usually be treated as elastic, at least 
for low electron energies. In addition, because the electron-neutral collisions are not long range, they can be 
assumed to be of very short duration (effectively instantaneous) compared to the time the electron spends between 
collisions. These assumptions allow one to derive an explicit form for the collision term. The derivation 
involves a consideration of the mechanics of electron scattering into and out of a specified region of velocity 
space. The appropriate collision term is shown below. 

Dcoii = n « J dv i J dn dn lv ' v 1 1 [f(v,) fn(vi) ' f(v) fn(vi)] 0> 

In Equation 7 f n is the distribution function of the background neutrals. n„ is the average density of these 

neutrals, vi defines the velocity space of the neutral particles and Q is a solid angle. ^ represents the 

differential cross section for the scattering processes under consideration and can be quite complex. The scattering 
cross section is generally a function of v'. Note the presence of the factor Iv-vjl in the integrand. This factor 
expresses the fact that no collisions can occur unless there is a relative velocity between the electrons and the 
background particles. The bracketed terms in the integrand represent two specific types of collisions. The first 
term, sometimes called a 'scattering in' term, defines those collisions which scatter electrons from all other 
velocities v' to the velocity of interest v. The second term, the 'scattering out' term, defines those collisions 
which scatter electrons from the velocity of interest v to any other velocity. Hence the 'scattering in' term is in a 
sense a source term in Equation 7, while the 'scattering out' term is a sink term. The balance of these terms 
along with the forcing terms from the left side of Equation 6 determines the behavior of f in a particular volume 
of velocity space. 

It can be shown that for a proper choice of f, the collision term can be made to vanish. This choice is a 
function only of v, and therefore in the absence of external forces (E, B, etc.) it is an equilibrium solution of the 
Boltzmann equation. This distribution is known as the Maxwell-Boltzmann distribution and is shown below. 

w.GfrJ’X-Ss) « 

In this expression T e is the electron temperature. It characterizes the size of the distribution in velocity space. 

Even given the explicit form for the collision term of Equation 7, the Boltzmann equation is very difficult 
to solve without further simplifications. In Section 2, a variety of simplifications to the collision term will be 
considered. In addition, there are some general assumptions that cover all work reported here. The Erst is that 
steady state will be assumed. The distribution function relaxes in time on the scale of several collision times. 
This will generally be on the order of a few tens of picoseconds for electrons in typical electric fields. 

A second assumption is that of homogeneity. That is, the electron distribution function will be assumed to 
be independent of x. This will be a reasonable approximation locally if the electric field is uniform over several 
electron mean free paths. The electron mean free path in air is approximately 10'6 m. Hence it takes at least 
1000 collision times for the electron to move a millimeter, even if it is assumed to move in a straight line. This 
time is much longer when the actual drift velocity of the electron distribution is used, rather than the straight line 
speed. Then if the electric field is uniform over fractions of millimeters the homogeneity approximation is 
acceptable. 
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The third basic assumption is that only electric field external forces are present. The effect of static 
magnetic fields that may be present will be ignored. The acceleration term under this assumption becomes 
eE 7 

a = - — ~ and Boltzmann’s equation reduces to, 
m 


§E JL = 3f\ 
m 3v z - dt/coll 


(9) 


Note here that E is assumed to be oriented along the positive z axis, but the z subscript will be suppressed. 


In Equation 9 f is now considered to be solely a function of v. This means that f can be written in general 
as fOx.Vy.Vj) or, given the cylindrical symmetry, f(v r ,v z ). The electric field E can be treated as a parameter, so 
there will be a family of distribution functions based on the value of E. 


The Maxwell-Boltzmann distribution given by Equation 8 represents an equilibrium solution of 
Boltzmann's equation in the absence of an electric field. This is die distribution to which collisions, and the 
collision term of Equation 9, drive the system. A simple way to think of the right hand side of the equation, 
therefore, is as a term tending to force the system back to equilibrium. Modeling this mathematically can be 
accomplished simply by replacing the full collision term with a relaxation process as shown below. 


3f\ f(») - fo 

dt/ coll “ * x 


( 10 ) 


Here f 0 is the Maxwell-Boltzmann equilibrium distribution, f(v) is the (perturbed) distribution in the presence of 
the electric field, and x is an average collision time for electron-neutral interactions, x can be replaced, if desired, 
by its inverse, the collision frequency v. The simplified Boltzmann equation then becomes, 

— ^- = v(f-Q (11) 

m 3v z w v 7 


The form of Equation 1 1 can be interpreted in much the same way as the more general Equation 7. The first term 
on the right hand side can be thought of as a term scattering electrons out of the volume dv centered on v. The 
second term represents scattering of electrons into the volume dv centered on v from all other velocities. 
Replacement of the full collision term with a relaxation term as in Equation 11 constitutes what is called the 
Krook model. 

2.1 LINEARIZED SOLUTION OF KROOK MODEL 

If the electric field strength is small enough, one expects the distribution function f to differ only slightly 
from the equilibrium distribution f 0 . To investigate this regime, f may be considered as the sum of f 0 and a 
(small) perturbation fy. 


f=f 0 + fi (12) 

This leads to a linearized solution for f in cylindrical coordinates. 

f(v„vj -fo + fl- (sjr) 3 " «P (• (l • f£) 03) 

Equation 13 is the solution to the linearized Krook model. Note in particular the last term in the equation. This 
term can be negative if E or v z are too large. Because f must be everywhere nonnegative, this represents a lowest 
order limit on the value of E which can be allowed in the linearized model. 

Figure 1 shows contours of the distribution function f as calculated from Equation 13 for E * 3 MV/m. 
The electron-neutral collision frequency, v, was chosen to be 1 x 10 12 s* 1 , and the electron temperature (of the 
unperturbed distribution f 0 ) was chosen to be 300* K. Thfe figure shows contours of f, with the dashed lines 
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indicating negative (nonphysical) values. The contours represent a projection of the full three dimensional 
distribution function. Rotation about the left side of the figure reproduces the third dimension. It should be 
noted that the linearized Krook model gives physically realistic solutions only for much lower values of electric 
field (~ 100 KV/m). 



V r (m/s x 10 5 ) 

Figure 1 Linearized Krook Model Distribution Function, E = 3 MV/m. Peak Value is 
1.11 x 10’ 1S 



Although linearization of Equation 11 leads to a simple solution in the Krook model, it is not necessary to 
linearize to obtain a solution. Equation 11 as it stands is a simple first order linear differential equation with 
constant coefficients and can be integrated directly. First, rewriting the equation in a more suitable form. 


9f mv . mv . 
3v z 'eE eE° 


(14) 


With the application of an integrating factor. Equation 14 becomes. 




In integrating the equation some care must be exercised in choosing limits of integration. If one integrates over 
v z from -oo to v z , the left hand side evaluated at may or may not vanish depending on the behavior of f at large 
negative v z . More convenient is to integrate from v z to +©°. At v z — > +°° the left hand side vanishes unless f 
grows faster exponentially than the integrating factor decays. But the distribution function f must go to zero as 
v z approaches infinity to be normalizable, so this is not a problem. Substituting for f 0 and performing the 
integration on the left hand side of Equation 15 leads to. 
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After some manipulation the closed form solution for fO^Vj) becomes. 


f(v r .vz) = 


m 2 v JT 
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It may be noted that Equation 17 is significantly more complex than Equation 13, the linearized Krook 
solution. Another difference is that f(v T ,v 2 ) as calculated from Equation 17 is everywhere positive, as required for 
distribution functions. 

Figure 2 is a contour plot of Equation 17 for electric field value of 3 MV/m. The figure is analogous to 
Figure 1 for the linearized solution. This solution is much more physically reasonable than that of Figure 1. 
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Figure 2 Krook Model Distribution Function, E = 3 MV/m 



2.3 DISCUSSION OF KROOK MODEL 

Because the Krook model is one of the easiest to deal with in the solution of the Boltzmann equation, it is 
of interest to briefly examine some implications of the approximations made. One of these is the form of the 
Boltzmann collision integral in the Krook model. In the Krook model t is defined to be a mean time between 
collisions, and can be related to a mean collision frequency, v. Defined as means, these quantities are not 
functions of velocity. Physically, however, they should be functions of velocity (or more accurately, speed), and 
t should be a function of velocity in Equation 10. The faster an electron moves through an essentially 
unmoving background, the more collisions it experiences in unit time. The mean free path of the electron may 
remain constant, but the time between collisions cannot. Then, because v should be monotonically increasing 
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with the speed of the electron, the relaxation term on the right of Equation 11 should be less important for 
electron speeds resulting in collision frequencies below the mean, and more important for electron speeds 
resulting in collision frequencies above the mean. This will tend to raise the distribution function f at low speeds 
and decrease it at higher speeds, leading to contours which appear more squashed than the extended distributions 
shown in Figure 2. 

A second issue has to do with normalization of the distribution. Recall from Equation 1 die normalization 
requirement on the function f, essentially that the integral of f over all velocity space must equal one. Consider 
now the Krook model of the Boltzmann equation (Equation 1 1), repeated below. 

First, assume that v is not a function of velocity, and integrate both sides of the equation over all v z (-°° to 
+oo). The left hand side is identically zero, because f is assumed to vanish as v z -> ±«. Passing v through the 
integrals on the right hand side leads to, 


J f dv z » J f 0 dv z (18) 

Then f will be properly normalized as long as f 0 is a normalized distribution, which it is when the Maxwellian 
distribution is used. 

Now consider what happens when v is allowed to be a function of electron speed. The left hand side still 
vanishes, for the same reason as before. The collision frequency on the right hand side, however, can no longer 
be passed through the integral, so one finds, 

Jv(lvl) f dv z * Jv(lvl) f 0 dv z (19) 

Now even though f 0 may be normalized. Equation 19 does not require the normalization of f. In general a 
normalized solution of Equation 11 cannot be found for this case, because the equation is nonlinear with the 
introduction of v(lvl). Hence the simple Krook model does not adequately handle the case of collision frequency 
which is a function of electron speed. 

A more realistic model would allow the collision frequency v to vary with the speed of the electron, v «= 
Ivl, or use the full Boltzmann collision integral. Both of these techniques were tried under this effort, but were 
found to lead to serious computational complexities beyond the scope of the effort. The interested reader is 
referred to Reference 4 for the details of these calculations. 

3.0 ELECTRON DISTRIBUTION FUNCTION CALCULATION USING 
MONTE CARLO TECHNIQUE 

The difficulty in dealing with any but the simplest cases in the Boltzmann collision term leads one to 
search for more convenient ways to calculate the electron distribution function. An obvious technique which has 
become practical only with the advent of digital computers is to simulate the motion of a large number of 
electrons. The simulation should include the effects of electric field force and interactions with the background 
neutral gas. Randomness is introduced into the simulation through the parameters of the collisions that the 
electrons undergo. 

The Monte Carlo simulation proceeds in the following way. An electron is dropped at t = 0 with no 
velocity into a medium containing an electric field. The electric field is assumed, without loss of generality, to 
be oriented along the positive z axis of a cylindrical coordinate system. The electron will be continuously 
accelerated by the electric field and will experience collisions with the background gas that randomly change its 
velocity. The distance the electron travels between collisions, 4 is based on its mean free path in the medium and 
is a random variable. Given 4 the next step in the simulation is to calculate the time At that the electron takes to 
travel that distance. At depends on the initial velocity of the election (its velocity immediately after the preceding 
collision), the distance 4 and the strength of the electric field. After the time At the electron is assumed to 
experience a collision with a background particle. The parameters of the collision, consisting of an impact 
parameter and an angle, are random variables. The collision calculation is the point at which a variety of 
physical models for electron/background interaction can be included. A simple model is considered here, in which 
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the electrons and background particles are assumed to behave as hard spheres. The sequence of steps above is 
continued for the given electron until a time to is reached, at which time the electron velocity is frozen in place. 
A second electron is then released with zero velocity and the sequence repeated. The time to is chosen so as to 
ensure that each electron has undergone a large number of collisions before its velocity is frozen. 

After a large number of electrons have been simulated (up to one million in the work reported here) a 
distribution function can be derived from the final calculated velocities of the electrons. Normalizations are 
necessary to account for the transformation from an intrinsically three dimensional electron motion to a two 
dimensional distribution function. 

The study of electron motions in Monte Carlo simulations has been considered by a number of researchers 
[5-8]. The technique has been used to study the growth of ionization in electron avalanches and electron 
diffusion. It is well respected in that far more realistic physics can be included than in more analytically oriented 
techniques. For example, it is possible in principle to include very complex variations in collision cross section 
with electron energy in a Monte Carlo calculation. This is extremely difficult in Boltzmann equation approaches. 
The present work is restricted to simple hard sphere collisions because of resource restrictions, not because more 
physics is difficult to include in the Monte Carlo framework. 

To demonstrate the Monte Carlo technique a simulation was undertaken to determine the electron 
distribution function in an airlike medium with an electric field strength of 300 kV/m. The value of 300 kV/m 
was chosen to avoid the physical process of ionization which in principle should be included for high intensity 
electric field cases. Ionization is of course a highly inelastic process and would invalidate the simple hard sphere 
scattering model used here. Inclusion of ionization processes can be accomplished through the use of additional 
random variables which define the probability of inelastic versus elastic collisions and the implementation of 
different scattering models for the two processes. In general the additional random variables need to account for 
varying scattering cross sections as a function of electron energy. 

As mentioned above, the simulation reported here used elastic hard sphere collisions only. A total of one 
million electrons were simulated. Each of the electrons was followed for 20 picoseconds (to), a time which 
allows on the average about twenty collisions with the background gas. Ideally, it would be desirable to have 
more collisions per electron. The simulation was truncated at 20 picoseconds to limit computer resource usage. 
The mean free path chosen for the simulation was 6.88 x 10' 7 m. At the end of the 20 picosecond simulation 
time the velocity of each electron was recorded and a new electron started. The distribution of final energies of all 
the electrons is shown in Figure 3. The distribution approximates quite well the analytically predicted 
distribution, known as a Druyvesteynian. The actual distribution function of the electrons is shown in Figure 4. 
The raggedness of the contours in Figure 4 requires some explanation. Although one million electrons were 
simulated, the vast majority of these electrons end up with final velocities having large r components. Required 
normalization reduces the importance of these electrons and increases the importance of those having small r 
components of velocity in the final distribution. This normalization, however, means that a relatively small 
number of electrons determines the electron velocity distribution function, resulting in statistical fluctuations in 
the contours of Figure 4. Adding more electrons to the simulation would eventually result in smoother contours, 
but an order of magnitude increase may be necessary. 

Comparison of the distribution of Figure 4 with that of Figure 2 reveals the same general behavior, as 
expected. More exact comparisons would require a larger number of electrons to be simulated, as well as a longer 
simulation time, both beyond the scope of this effort. For simple cases in which a Boltzmann solution is 
possible, it is clear that the Boltzmann technique is both faster and simpler. The main advantage of the Monte 
Carlo technique is in its ability to add complex physical processes without correspondingly complex increases in 
difficulty. That is, the Monte Carlo technique takes significant computer resources to generate a distribution 
function even in simple cases, but the extension to more realistic physics is not difficult or significantly more 
costly. Solutions of the Boltzmann equation are fast and inexpensive for simple cases, but very difficult for more 
realistic cases. 
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Figure 3 Distribution of Electron Energies at Conclusion of Monte Carlo Simulation of 
One Million Electrons. Tbe Fundamental Form of the Distribution is Known as 
a Druyvestinian 



Figure 4 Electron Distribution Function From Monte Carlo Simulation of One Million 
Electrons. The Function is Unnormalized 
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